// Copyright (C) 2008, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#include "CoinUtilsConfig.h"
#include "CoinHelperFunctions.hpp"
#include "CoinStructuredModel.hpp"
#include "CoinSort.hpp"
#include "CoinMpsIO.hpp"
#include "CoinMessage.hpp"
#include "CoinFloatEqual.hpp"

//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
CoinStructuredModel::CoinStructuredModel()
  : CoinBaseModel()
  , numberRowBlocks_(0)
  , numberColumnBlocks_(0)
  , numberElementBlocks_(0)
  , maximumElementBlocks_(0)
  , blocks_(NULL)
  , coinModelBlocks_(NULL)
  , blockType_(NULL)
{
}
/* Read a problem in MPS or GAMS format from the given filename.
 */
CoinStructuredModel::CoinStructuredModel(const char *fileName,
  int decomposeType,
  int maxBlocks)
  : CoinBaseModel()
  , numberRowBlocks_(0)
  , numberColumnBlocks_(0)
  , numberElementBlocks_(0)
  , maximumElementBlocks_(0)
  , blocks_(NULL)
  , coinModelBlocks_(NULL)
  , blockType_(NULL)
{
  CoinModel coinModel(fileName, false);
  if (coinModel.numberRows()) {
    problemName_ = coinModel.getProblemName();
    optimizationDirection_ = coinModel.optimizationDirection();
    objectiveOffset_ = coinModel.objectiveOffset();
    if (!decomposeType) {
      addBlock("row_master", "column_master", coinModel);
    } else {
      const CoinPackedMatrix *matrix = coinModel.packedMatrix();
      if (!matrix)
        coinModel.convertMatrix();
      decompose(coinModel, decomposeType, maxBlocks);
      //addBlock("row_master","column_master",coinModel);
    }
  }
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
CoinStructuredModel::CoinStructuredModel(const CoinStructuredModel &rhs)
  : CoinBaseModel(rhs)
  , numberRowBlocks_(rhs.numberRowBlocks_)
  , numberColumnBlocks_(rhs.numberColumnBlocks_)
  , numberElementBlocks_(rhs.numberElementBlocks_)
  , maximumElementBlocks_(rhs.maximumElementBlocks_)
{
  if (maximumElementBlocks_) {
    blocks_ = CoinCopyOfArray(rhs.blocks_, maximumElementBlocks_);
    for (int i = 0; i < numberElementBlocks_; i++)
      blocks_[i] = rhs.blocks_[i]->clone();
    blockType_ = CoinCopyOfArray(rhs.blockType_, maximumElementBlocks_);
    if (rhs.coinModelBlocks_) {
      coinModelBlocks_ = CoinCopyOfArray(rhs.coinModelBlocks_,
        maximumElementBlocks_);
      for (int i = 0; i < numberElementBlocks_; i++)
        coinModelBlocks_[i] = new CoinModel(*rhs.coinModelBlocks_[i]);
    } else {
      coinModelBlocks_ = NULL;
    }
  } else {
    blocks_ = NULL;
    blockType_ = NULL;
    coinModelBlocks_ = NULL;
  }
  rowBlockNames_ = rhs.rowBlockNames_;
  columnBlockNames_ = rhs.columnBlockNames_;
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
CoinStructuredModel::~CoinStructuredModel()
{
  for (int i = 0; i < numberElementBlocks_; i++)
    delete blocks_[i];
  delete[] blocks_;
  delete[] blockType_;
  if (coinModelBlocks_) {
    for (int i = 0; i < numberElementBlocks_; i++)
      delete coinModelBlocks_[i];
    delete[] coinModelBlocks_;
  }
}

// Clone
CoinBaseModel *
CoinStructuredModel::clone() const
{
  return new CoinStructuredModel(*this);
}

//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
CoinStructuredModel &
CoinStructuredModel::operator=(const CoinStructuredModel &rhs)
{
  if (this != &rhs) {
    CoinBaseModel::operator=(rhs);
    for (int i = 0; i < numberElementBlocks_; i++)
      delete blocks_[i];
    delete[] blocks_;
    delete[] blockType_;
    if (coinModelBlocks_) {
      for (int i = 0; i < numberElementBlocks_; i++)
        delete coinModelBlocks_[i];
      delete[] coinModelBlocks_;
    }
    numberRowBlocks_ = rhs.numberRowBlocks_;
    numberColumnBlocks_ = rhs.numberColumnBlocks_;
    numberElementBlocks_ = rhs.numberElementBlocks_;
    maximumElementBlocks_ = rhs.maximumElementBlocks_;
    if (maximumElementBlocks_) {
      blocks_ = CoinCopyOfArray(rhs.blocks_, maximumElementBlocks_);
      for (int i = 0; i < numberElementBlocks_; i++)
        blocks_[i] = rhs.blocks_[i]->clone();
      blockType_ = CoinCopyOfArray(rhs.blockType_, maximumElementBlocks_);
      if (rhs.coinModelBlocks_) {
        coinModelBlocks_ = CoinCopyOfArray(rhs.coinModelBlocks_,
          maximumElementBlocks_);
        for (int i = 0; i < numberElementBlocks_; i++)
          coinModelBlocks_[i] = new CoinModel(*rhs.coinModelBlocks_[i]);
      } else {
        coinModelBlocks_ = NULL;
      }
    } else {
      blocks_ = NULL;
      blockType_ = NULL;
      coinModelBlocks_ = NULL;
    }
    rowBlockNames_ = rhs.rowBlockNames_;
    columnBlockNames_ = rhs.columnBlockNames_;
  }
  return *this;
}
static bool sameValues(const double *a, const double *b, int n)
{
  int i;
  for (i = 0; i < n; i++) {
    if (a[i] != b[i])
      break;
  }
  return (i == n);
}
static bool sameValues(const int *a, const int *b, int n)
{
  int i;
  for (i = 0; i < n; i++) {
    if (a[i] != b[i])
      break;
  }
  return (i == n);
}
static bool sameValues(const CoinModel *a, const CoinModel *b, bool doRows)
{
  int i = 0;
  int n = 0;
  if (doRows) {
    n = a->numberRows();
    for (i = 0; i < n; i++) {
      const char *aName = a->getRowName(i);
      const char *bName = b->getRowName(i);
      bool good = true;
      if (aName) {
        if (!bName || strcmp(aName, bName))
          good = false;
      } else if (bName) {
        good = false;
      }
      if (!good)
        break;
    }
  } else {
    n = a->numberColumns();
    for (i = 0; i < n; i++) {
      const char *aName = a->getColumnName(i);
      const char *bName = b->getColumnName(i);
      bool good = true;
      if (aName) {
        if (!bName || strcmp(aName, bName))
          good = false;
      } else if (bName) {
        good = false;
      }
      if (!good)
        break;
    }
  }
  return (i == n);
}
// Add a row block name and number of rows
int CoinStructuredModel::addRowBlock(int numberRows, const std::string &name)
{
  int iRowBlock;
  for (iRowBlock = 0; iRowBlock < numberRowBlocks_; iRowBlock++) {
    if (name == rowBlockNames_[iRowBlock])
      break;
  }
  if (iRowBlock == numberRowBlocks_) {
    rowBlockNames_.push_back(name);
    numberRowBlocks_++;
    numberRows_ += numberRows;
  }
  return iRowBlock;
}
// Return a row block index given a row block name
int CoinStructuredModel::rowBlock(const std::string &name) const
{
  int iRowBlock;
  for (iRowBlock = 0; iRowBlock < numberRowBlocks_; iRowBlock++) {
    if (name == rowBlockNames_[iRowBlock])
      break;
  }
  if (iRowBlock == numberRowBlocks_)
    iRowBlock = -1;
  return iRowBlock;
}
// Add a column block name and number of columns
int CoinStructuredModel::addColumnBlock(int numberColumns, const std::string &name)
{
  int iColumnBlock;
  for (iColumnBlock = 0; iColumnBlock < numberColumnBlocks_; iColumnBlock++) {
    if (name == columnBlockNames_[iColumnBlock])
      break;
  }
  if (iColumnBlock == numberColumnBlocks_) {
    columnBlockNames_.push_back(name);
    numberColumnBlocks_++;
    numberColumns_ += numberColumns;
  }
  return iColumnBlock;
}
// Return a column block index given a column block name
int CoinStructuredModel::columnBlock(const std::string &name) const
{
  int iColumnBlock;
  for (iColumnBlock = 0; iColumnBlock < numberColumnBlocks_; iColumnBlock++) {
    if (name == columnBlockNames_[iColumnBlock])
      break;
  }
  if (iColumnBlock == numberColumnBlocks_)
    iColumnBlock = -1;
  return iColumnBlock;
}
// Return number of elements
CoinBigIndex
CoinStructuredModel::numberElements() const
{
  CoinBigIndex numberElements = 0;
  for (int iBlock = 0; iBlock < numberElementBlocks_; iBlock++) {
    numberElements += blocks_[iBlock]->numberElements();
  }
  return numberElements;
}
// Return i'th block as CoinModel (or NULL)
CoinModel *
CoinStructuredModel::coinBlock(int i) const
{
  CoinModel *block = dynamic_cast< CoinModel * >(blocks_[i]);
  if (block)
    return block;
  else if (coinModelBlocks_)
    return coinModelBlocks_[i];
  else
    return NULL;
}

/* Return block as a CoinModel block
   and fill in info structure and update counts
*/
CoinModel *
CoinStructuredModel::coinModelBlock(CoinModelBlockInfo &info)
{
  // CoinStructuredModel
  int numberRows = this->numberRows();
  int numberColumns = this->numberColumns();
  int numberRowBlocks = this->numberRowBlocks();
  int numberColumnBlocks = this->numberColumnBlocks();
  int numberElementBlocks = this->numberElementBlocks();
  CoinBigIndex numberElements = this->numberElements();
  // See what is needed
  double *rowLower = NULL;
  double *rowUpper = NULL;
  double *columnLower = NULL;
  double *columnUpper = NULL;
  double *objective = NULL;
  int *integerType = NULL;
  info = CoinModelBlockInfo();
  CoinModel **blocks = new CoinModel *[numberElementBlocks];
  for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
    CoinModelBlockInfo thisInfo = blockType_[iBlock];
    CoinStructuredModel *subModel = dynamic_cast< CoinStructuredModel * >(blocks_[iBlock]);
    CoinModel *thisBlock;
    if (subModel) {
      thisBlock = subModel->coinModelBlock(thisInfo);
      fillInfo(thisInfo, subModel);
      setCoinModel(thisBlock, iBlock);
    } else {
      thisBlock = dynamic_cast< CoinModel * >(blocks_[iBlock]);
      assert(thisBlock);
      fillInfo(thisInfo, thisBlock);
    }
    blocks[iBlock] = thisBlock;
    if (thisInfo.rhs && !info.rhs) {
      info.rhs = 1;
      rowLower = new double[numberRows];
      rowUpper = new double[numberRows];
      CoinFillN(rowLower, numberRows, -COIN_DBL_MAX);
      CoinFillN(rowUpper, numberRows, COIN_DBL_MAX);
    }
    if (thisInfo.bounds && !info.bounds) {
      info.bounds = 1;
      columnLower = new double[numberColumns];
      columnUpper = new double[numberColumns];
      objective = new double[numberColumns];
      CoinFillN(columnLower, numberColumns, 0.0);
      CoinFillN(columnUpper, numberColumns, COIN_DBL_MAX);
      CoinFillN(objective, numberColumns, 0.0);
    }
    if (thisInfo.integer && !info.integer) {
      info.integer = 1;
      integerType = new int[numberColumns];
      CoinFillN(integerType, numberColumns, 0);
    }
    if (thisInfo.rowName && !info.rowName) {
      info.rowName = 1;
    }
    if (thisInfo.columnName && !info.columnName) {
      info.columnName = 1;
    }
  }
  // Space for elements
  int *row = new int[numberElements];
  int *column = new int[numberElements];
  double *element = new double[numberElements];
  numberElements = 0;
  // Bases for blocks
  int *rowBase = new int[numberRowBlocks];
  CoinFillN(rowBase, numberRowBlocks, -1);
  CoinModelBlockInfo *rowBlockInfo = new CoinModelBlockInfo[numberRowBlocks];
  int *columnBase = new int[numberColumnBlocks];
  CoinFillN(columnBase, numberColumnBlocks, -1);
  CoinModelBlockInfo *columnBlockInfo = new CoinModelBlockInfo[numberColumnBlocks];
  for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
    int iRowBlock = rowBlock(blocks[iBlock]->getRowBlock());
    assert(iRowBlock >= 0 && iRowBlock < numberRowBlocks);
    if (rowBase[iRowBlock] == -1)
      rowBase[iRowBlock] = blocks[iBlock]->numberRows();
    else
      assert(rowBase[iRowBlock] == blocks[iBlock]->numberRows());
    int iColumnBlock = columnBlock(blocks[iBlock]->getColumnBlock());
    assert(iColumnBlock >= 0 && iColumnBlock < numberColumnBlocks);
    if (columnBase[iColumnBlock] == -1)
      columnBase[iColumnBlock] = blocks[iBlock]->numberColumns();
    else
      assert(columnBase[iColumnBlock] == blocks[iBlock]->numberColumns());
  }
  int n = 0;
  for (int iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
    int k = rowBase[iBlock];
    rowBase[iBlock] = n;
    assert(k >= 0);
    n += k;
  }
  assert(n == numberRows);
  n = 0;
  for (int iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
    int k = columnBase[iBlock];
    columnBase[iBlock] = n;
    assert(k >= 0);
    n += k;
  }
  assert(n == numberColumns);
  for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
    CoinModelBlockInfo thisInfo = blockType_[iBlock];
    CoinModel *thisBlock = blocks[iBlock];
    int iRowBlock = rowBlock(blocks[iBlock]->getRowBlock());
    int iRowBase = rowBase[iRowBlock];
    int nRows = thisBlock->numberRows();
    // could check if identical before error
    if (thisInfo.rhs) {
      assert(!rowBlockInfo[iRowBlock].rhs);
      rowBlockInfo[iRowBlock].rhs = 1;
      memcpy(rowLower + iRowBase, thisBlock->rowLowerArray(),
        nRows * sizeof(double));
      memcpy(rowUpper + iRowBase, thisBlock->rowUpperArray(),
        nRows * sizeof(double));
    }
    int iColumnBlock = columnBlock(blocks[iBlock]->getColumnBlock());
    int iColumnBase = columnBase[iColumnBlock];
    int nColumns = thisBlock->numberColumns();
    if (thisInfo.bounds) {
      assert(!columnBlockInfo[iColumnBlock].bounds);
      columnBlockInfo[iColumnBlock].bounds = 1;
      memcpy(columnLower + iColumnBase, thisBlock->columnLowerArray(),
        nColumns * sizeof(double));
      memcpy(columnUpper + iColumnBase, thisBlock->columnUpperArray(),
        nColumns * sizeof(double));
      memcpy(objective + iColumnBase, thisBlock->objectiveArray(),
        nColumns * sizeof(double));
    }
    if (thisInfo.integer) {
      assert(!columnBlockInfo[iColumnBlock].integer);
      columnBlockInfo[iColumnBlock].integer = 1;
      memcpy(integerType + iColumnBase, thisBlock->integerTypeArray(),
        nColumns * sizeof(int));
    }
    const CoinPackedMatrix *elementBlock = thisBlock->packedMatrix();
    // get matrix data pointers
    const int *row2 = elementBlock->getIndices();
    const CoinBigIndex *columnStart = elementBlock->getVectorStarts();
    const double *elementByColumn = elementBlock->getElements();
    const int *columnLength = elementBlock->getVectorLengths();
    int n = elementBlock->getNumCols();
    assert(elementBlock->isColOrdered());
    for (int iColumn = 0; iColumn < n; iColumn++) {
      CoinBigIndex j;
      int jColumn = iColumn + iColumnBase;
      for (j = columnStart[iColumn];
           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
        row[numberElements] = row2[j] + iRowBase;
        column[numberElements] = jColumn;
        element[numberElements++] = elementByColumn[j];
      }
    }
  }
  delete[] rowBlockInfo;
  delete[] columnBlockInfo;
  CoinPackedMatrix matrix(true, row, column, element, numberElements);
  if (numberElements)
    info.matrix = 1;
  delete[] row;
  delete[] column;
  delete[] element;
  CoinModel *block = new CoinModel(numberRows, numberColumns,
    &matrix,
    rowLower, rowUpper,
    columnLower, columnUpper, objective);
  delete[] rowLower;
  delete[] rowUpper;
  delete[] columnLower;
  delete[] columnUpper;
  delete[] objective;
  // Do integers if wanted
  if (integerType) {
    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
      block->setColumnIsInteger(iColumn, integerType[iColumn] != 0);
    }
  }
  delete[] integerType;
  block->setObjectiveOffset(objectiveOffset());
  if (info.rowName || info.columnName) {
    for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
      CoinModelBlockInfo thisInfo;
      CoinModel *thisBlock = blocks[iBlock];
      int iRowBlock = rowBlock(thisBlock->getRowBlock());
      int iRowBase = rowBase[iRowBlock];
      if (thisInfo.rowName) {
        int numberItems = thisBlock->rowNames()->numberItems();
        assert(thisBlock->numberRows() >= numberItems);
        if (numberItems) {
          const char *const *rowNames = thisBlock->rowNames()->names();
          for (int i = 0; i < numberItems; i++) {
            std::string name = rowNames[i];
            block->setRowName(i + iRowBase, name.c_str());
          }
        }
      }
      int iColumnBlock = columnBlock(thisBlock->getColumnBlock());
      int iColumnBase = columnBase[iColumnBlock];
      if (thisInfo.columnName) {
        int numberItems = thisBlock->columnNames()->numberItems();
        assert(thisBlock->numberColumns() >= numberItems);
        if (numberItems) {
          const char *const *columnNames = thisBlock->columnNames()->names();
          for (int i = 0; i < numberItems; i++) {
            std::string name = columnNames[i];
            block->setColumnName(i + iColumnBase, name.c_str());
          }
        }
      }
    }
  }
  delete[] rowBase;
  delete[] columnBase;
  for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
    if (static_cast< CoinBaseModel * >(blocks[iBlock]) != static_cast< CoinBaseModel * >(blocks_[iBlock]))
      delete blocks[iBlock];
  }
  delete[] blocks;
  return block;
}
// Sets given block into coinModelBlocks_
void CoinStructuredModel::setCoinModel(CoinModel *block, int iBlock)
{
  if (!coinModelBlocks_) {
    coinModelBlocks_ = new CoinModel *[maximumElementBlocks_];
    CoinZeroN(coinModelBlocks_, maximumElementBlocks_);
  }
  delete coinModelBlocks_[iBlock];
  coinModelBlocks_[iBlock] = block;
}
// Refresh info in blockType_
void CoinStructuredModel::refresh(int iBlock)
{
  fillInfo(blockType_[iBlock], coinBlock(iBlock));
}
/* Fill in info structure and update counts
   Returns number of inconsistencies on border
*/
int CoinStructuredModel::fillInfo(CoinModelBlockInfo &info,
  const CoinModel *block)
{
  int whatsSet = block->whatIsSet();
  info.matrix = static_cast< char >(((whatsSet & 1) != 0) ? 1 : 0);
  info.rhs = static_cast< char >(((whatsSet & 2) != 0) ? 1 : 0);
  info.rowName = static_cast< char >(((whatsSet & 4) != 0) ? 1 : 0);
  info.integer = static_cast< char >(((whatsSet & 32) != 0) ? 1 : 0);
  info.bounds = static_cast< char >(((whatsSet & 8) != 0) ? 1 : 0);
  info.columnName = static_cast< char >(((whatsSet & 16) != 0) ? 1 : 0);
  int numberRows = block->numberRows();
  int numberColumns = block->numberColumns();
  // Which block
  int iRowBlock = addRowBlock(numberRows, block->getRowBlock());
  info.rowBlock = iRowBlock;
  int iColumnBlock = addColumnBlock(numberColumns, block->getColumnBlock());
  info.columnBlock = iColumnBlock;
  int numberErrors = 0;
  CoinModelBlockInfo sumInfo = blockType_[numberElementBlocks_ - 1];
  int iRhs = (sumInfo.rhs) ? numberElementBlocks_ - 1 : -1;
  int iRowName = (sumInfo.rowName) ? numberElementBlocks_ - 1 : -1;
  int iBounds = (sumInfo.bounds) ? numberElementBlocks_ - 1 : -1;
  int iColumnName = (sumInfo.columnName) ? numberElementBlocks_ - 1 : -1;
  int iInteger = (sumInfo.integer) ? numberElementBlocks_ - 1 : -1;
  for (int i = 0; i < numberElementBlocks_ - 1; i++) {
    if (iRowBlock == blockType_[i].rowBlock) {
      if (numberRows != blocks_[i]->numberRows())
        numberErrors += 1000;
      if (blockType_[i].rhs) {
        if (iRhs < 0) {
          iRhs = i;
        } else {
          // check
          const double *a = static_cast< CoinModel * >(blocks_[iRhs])->rowLowerArray();
          const double *b = static_cast< CoinModel * >(blocks_[i])->rowLowerArray();
          if (!sameValues(a, b, numberRows))
            numberErrors++;
          a = static_cast< CoinModel * >(blocks_[iRhs])->rowUpperArray();
          b = static_cast< CoinModel * >(blocks_[i])->rowUpperArray();
          if (!sameValues(a, b, numberRows))
            numberErrors++;
        }
      }
      if (blockType_[i].rowName) {
        if (iRowName < 0) {
          iRowName = i;
        } else {
          // check
          if (!sameValues(static_cast< CoinModel * >(blocks_[iRowName]),
                static_cast< CoinModel * >(blocks_[i]), true))
            numberErrors++;
        }
      }
    }
    if (iColumnBlock == blockType_[i].columnBlock) {
      if (numberColumns != blocks_[i]->numberColumns())
        numberErrors += 1000;
      if (blockType_[i].bounds) {
        if (iBounds < 0) {
          iBounds = i;
        } else {
          // check
          const double *a = static_cast< CoinModel * >(blocks_[iBounds])->columnLowerArray();
          const double *b = static_cast< CoinModel * >(blocks_[i])->columnLowerArray();
          if (!sameValues(a, b, numberColumns))
            numberErrors++;
          a = static_cast< CoinModel * >(blocks_[iBounds])->columnUpperArray();
          b = static_cast< CoinModel * >(blocks_[i])->columnUpperArray();
          if (!sameValues(a, b, numberColumns))
            numberErrors++;
          a = static_cast< CoinModel * >(blocks_[iBounds])->objectiveArray();
          b = static_cast< CoinModel * >(blocks_[i])->objectiveArray();
          if (!sameValues(a, b, numberColumns))
            numberErrors++;
        }
      }
      if (blockType_[i].columnName) {
        if (iColumnName < 0) {
          iColumnName = i;
        } else {
          // check
          if (!sameValues(static_cast< CoinModel * >(blocks_[iColumnName]),
                static_cast< CoinModel * >(blocks_[i]), false))
            numberErrors++;
        }
      }
      if (blockType_[i].integer) {
        if (iInteger < 0) {
          iInteger = i;
        } else {
          // check
          const int *a = static_cast< CoinModel * >(blocks_[iInteger])->integerTypeArray();
          const int *b = static_cast< CoinModel * >(blocks_[i])->integerTypeArray();
          if (!sameValues(a, b, numberColumns))
            numberErrors++;
        }
      }
    }
  }
  return numberErrors;
}
/* Fill in info structure and update counts
*/
void CoinStructuredModel::fillInfo(CoinModelBlockInfo &info,
  const CoinStructuredModel *block)
{
  int numberRows = block->numberRows();
  int numberColumns = block->numberColumns();
  // Which block
  int iRowBlock = addRowBlock(numberRows, block->rowBlockName_);
  info.rowBlock = iRowBlock;
  int iColumnBlock = addColumnBlock(numberColumns, block->columnBlockName_);
  info.columnBlock = iColumnBlock;
}
/* add a block from a CoinModel without names*/
int CoinStructuredModel::addBlock(const std::string &rowBlock,
  const std::string &columnBlock,
  const CoinBaseModel &block)
{
  CoinBaseModel *block2 = block.clone();
  return addBlock(rowBlock, columnBlock, block2);
}
/* add a block using names 
 */
int CoinStructuredModel::addBlock(const std::string &rowBlock,
  const std::string &columnBlock,
  const CoinPackedMatrix &matrix,
  const double *rowLower, const double *rowUpper,
  const double *columnLower, const double *columnUpper,
  const double *objective)
{
  CoinModel *block = new CoinModel();
  block->loadBlock(matrix, columnLower, columnUpper, objective,
    rowLower, rowUpper);
  return addBlock(rowBlock, columnBlock, block);
}
/* add a block from a CoinModel without names*/
int CoinStructuredModel::addBlock(const std::string &rowBlock,
  const std::string &columnBlock,
  CoinBaseModel *block)
{
  if (numberElementBlocks_ == maximumElementBlocks_) {
    maximumElementBlocks_ = 3 * (maximumElementBlocks_ + 10) / 2;
    CoinBaseModel **temp = new CoinBaseModel *[maximumElementBlocks_];
    memcpy(temp, blocks_, numberElementBlocks_ * sizeof(CoinBaseModel *));
    delete[] blocks_;
    blocks_ = temp;
    CoinModelBlockInfo *temp2 = new CoinModelBlockInfo[maximumElementBlocks_];
    memcpy(temp2, blockType_, numberElementBlocks_ * sizeof(CoinModelBlockInfo));
    delete[] blockType_;
    blockType_ = temp2;
    if (coinModelBlocks_) {
      CoinModel **temp = new CoinModel *[maximumElementBlocks_];
      CoinZeroN(temp, maximumElementBlocks_);
      memcpy(temp, coinModelBlocks_, numberElementBlocks_ * sizeof(CoinModel *));
      delete[] coinModelBlocks_;
      coinModelBlocks_ = temp;
    }
  }
  blocks_[numberElementBlocks_++] = block;
  block->setRowBlock(rowBlock);
  block->setColumnBlock(columnBlock);
  int numberErrors = 0;
  CoinModel *coinBlock = dynamic_cast< CoinModel * >(block);
  if (coinBlock) {
    // Convert matrix
    if (coinBlock->type() != 3)
      coinBlock->convertMatrix();
    numberErrors = fillInfo(blockType_[numberElementBlocks_ - 1], coinBlock);
  } else {
    CoinStructuredModel *subModel = dynamic_cast< CoinStructuredModel * >(block);
    assert(subModel);
    CoinModel *blockX = subModel->coinModelBlock(blockType_[numberElementBlocks_ - 1]);
    fillInfo(blockType_[numberElementBlocks_ - 1], subModel);
    setCoinModel(blockX, numberElementBlocks_ - 1);
  }
  return numberErrors;
}

/* add a block from a CoinModel with names*/
int CoinStructuredModel::addBlock(const CoinBaseModel &block)
{

  //inline const std::string & getRowBlock() const
  //abort();
  return addBlock(block.getRowBlock(), block.getColumnBlock(),
    block);
}
/* Decompose a model specified as arrays + CoinPackedMatrix
   1 - try D-W
   2 - try Benders
   3 - try Staircase
   Returns number of blocks or zero if no structure
*/
int CoinStructuredModel::decompose(const CoinPackedMatrix &matrix,
  const double *rowLower, const double *rowUpper,
  const double *columnLower, const double *columnUpper,
  const double *objective, int type, int maxBlocks,
  int *starts, double objectiveOffset)
{
  setObjectiveOffset(objectiveOffset);
  int numberBlocks = 0;
  char generalPrint[200];
  bool wantDecomposition = type > 2;
  type %= 10;
  if (type == 1) { // Try master at top and bottom
    // get row copy
    CoinPackedMatrix rowCopy = matrix;
    rowCopy.reverseOrdering();
    const int *row = matrix.getIndices();
    const int *columnLength = matrix.getVectorLengths();
    const CoinBigIndex *columnStart = matrix.getVectorStarts();
    //const double * elementByColumn = matrix.getElements();
    const int *column = rowCopy.getIndices();
    const int *rowLength = rowCopy.getVectorLengths();
    const CoinBigIndex *rowStart = rowCopy.getVectorStarts();
    //const double * elementByRow = rowCopy.getElements();
    int numberRows = matrix.getNumRows();
    int *rowBlock = new int[numberRows + 1];
    int iRow;
    // Row counts (maybe look at long rows for master)
    CoinZeroN(rowBlock, numberRows + 1);
    for (iRow = 0; iRow < numberRows; iRow++) {
      int length = rowLength[iRow];
      rowBlock[length]++;
    }
    if (!starts) {
      for (iRow = 0; iRow < numberRows + 1; iRow++) {
        if (rowBlock[iRow]) {
          sprintf(generalPrint, "%d rows have %d elements", rowBlock[iRow], iRow);
          handler_->message(COIN_GENERAL_INFO2, messages_) << generalPrint << CoinMessageEol;
        }
      }
    }
    bool newWay = true;
    // to say if column looked at
    int numberColumns = matrix.getNumCols();
    int *columnBlock = new int[numberColumns];
    int iColumn;
    int *whichRow = new int[numberRows];
    int *whichColumn = new int[numberColumns];
    int *stack = new int[numberRows];
    if (newWay && !starts) {
      //double best2[3]={COIN_DBL_MAX,COIN_DBL_MAX,COIN_DBL_MAX};
      double best2[3] = { 0.0, 0.0, 0.0 };
      int row2[3] = { -1, -1, -1 };
      // try forward and backward and sorted
      for (int iWay = 0; iWay < 3; iWay++) {
        if (iWay == 0) {
          // forwards
          for (int i = 0; i < numberRows; i++)
            stack[i] = i;
        } else if (iWay == 1) {
          // backwards
          for (int i = 0; i < numberRows; i++)
            stack[i] = numberRows - 1 - i;
        } else {
          // sparsest first
          for (int i = 0; i < numberRows; i++) {
            rowBlock[i] = rowLength[i];
            stack[i] = i;
          }
          CoinSort_2(rowBlock, rowBlock + numberRows, stack);
        }
        CoinFillN(rowBlock, numberRows, -1);
        rowBlock[numberRows] = 0;
        CoinFillN(whichRow, numberRows, -1);
        CoinFillN(columnBlock, numberColumns, -1);
        CoinFillN(whichColumn, numberColumns, -1);
        int numberMarkedColumns = 0;
        numberBlocks = 0;
        int bestRow = -1;
        int maximumInBlock = 0;
        int rowsDone = 0;
        int checkAfterRows = (5 * numberRows) / 10 + 1;
#define OSL_WAY
#ifdef OSL_WAY
        double best = COIN_DBL_MAX;
        int bestRowsDone = -1;
#else
        double best = 0.0; //COIN_DBL_MAX;
#endif
        int numberGoodBlocks = 0;

        for (int kRow = 0; kRow < numberRows; kRow++) {
          iRow = stack[kRow];
          CoinBigIndex start = rowStart[iRow];
          CoinBigIndex end = start + rowLength[iRow];
          int iBlock = -1;
          for (CoinBigIndex j = start; j < end; j++) {
            int iColumn = column[j];
            if (columnBlock[iColumn] >= 0) {
              // already marked
              if (iBlock < 0) {
                iBlock = columnBlock[iColumn];
              } else if (iBlock != columnBlock[iColumn]) {
                // join two blocks
                int jBlock = columnBlock[iColumn];
                numberGoodBlocks--;
                // Increase count of iBlock
                rowBlock[iBlock] += rowBlock[jBlock];
                rowBlock[jBlock] = 0;
                // First column of block jBlock
                int jColumn = whichRow[jBlock];
                while (jColumn >= 0) {
                  columnBlock[jColumn] = iBlock;
                  iColumn = jColumn;
                  jColumn = whichColumn[jColumn];
                }
                whichColumn[iColumn] = whichRow[iBlock];
                whichRow[iBlock] = whichRow[jBlock];
                whichRow[jBlock] = -1;
              }
            }
          }
          int n = static_cast< int >(end - start);
          // If not in block - then start one
          if (iBlock < 0) {
            // unless null row
            if (n) {
              iBlock = numberBlocks;
              numberBlocks++;
              numberGoodBlocks++;
              int jColumn = column[start];
              columnBlock[jColumn] = iBlock;
              whichRow[iBlock] = jColumn;
              numberMarkedColumns += n;
              rowBlock[iBlock] = n;
              for (CoinBigIndex j = start + 1; j < end; j++) {
                int iColumn = column[j];
                columnBlock[iColumn] = iBlock;
                whichColumn[jColumn] = iColumn;
                jColumn = iColumn;
              }
            } else {
              rowBlock[numberRows]++;
            }
          } else {
            // add all to this block if not already in
            int jColumn = whichRow[iBlock];
            for (CoinBigIndex j = start; j < end; j++) {
              int iColumn = column[j];
              if (columnBlock[iColumn] < 0) {
                numberMarkedColumns++;
                rowBlock[iBlock]++;
                columnBlock[iColumn] = iBlock;
                whichColumn[iColumn] = jColumn;
                jColumn = iColumn;
              }
            }
            whichRow[iBlock] = jColumn;
          }
#if 0
	  {
	    int nn=0;
	    int * temp = new int [numberRows];
	    CoinZeroN(temp,numberRows);
	    for (int i=0;i<numberColumns;i++) {
	      int iBlock = columnBlock[i];
	      if (iBlock>=0) {
		nn++;
		assert (iBlock<numberBlocks);
		temp[iBlock]++;
	      }
	    }
	    for (int i=0;i<numberBlocks;i++)
	      assert (temp[i]==rowBlock[i]);
	    assert (nn==numberMarkedColumns);
	    for (int i=0;i<numberBlocks;i++) {
	      // First column of block i
	      int jColumn = whichRow[i];
	      int n=0;
	      while (jColumn>=0) {
		n++;
		jColumn = whichColumn[jColumn];
	      }
	      assert (n==temp[i]);
	    }
	    delete [] temp;
	  }
#endif
          rowsDone++;
          if (iBlock >= 0)
            maximumInBlock = std::max(maximumInBlock, rowBlock[iBlock]);
          if (rowsDone >= checkAfterRows) {
            assert(numberGoodBlocks > 0);
            double averageSize = static_cast< double >(numberMarkedColumns) / static_cast< double >(numberGoodBlocks);
#ifndef OSL_WAY
            double notionalBlocks = static_cast< double >(numberMarkedColumns) / averageSize;
            if (maximumInBlock < 3 * averageSize && numberGoodBlocks > 2) {
              if (best * (numberRows - rowsDone) < notionalBlocks) {
                best = notionalBlocks / static_cast< double >(numberRows - rowsDone);
                bestRow = kRow;
              }
            }
#else
            if (maximumInBlock * 10 < numberColumns * 11 && numberGoodBlocks > 1) {
              double test = maximumInBlock + 0.0 * averageSize;
              if (best * static_cast< double >(rowsDone) > test) {
                best = test / static_cast< double >(rowsDone);
                bestRow = kRow;
                bestRowsDone = rowsDone;
              }
            }
#endif
          }
        }
#ifndef OSL_WAY
        best2[iWay] = best;
#else
        if (bestRowsDone < numberRows)
          best2[iWay] = -(numberRows - bestRowsDone);
        else
          best2[iWay] = -numberRows;
#endif
        row2[iWay] = bestRow;
      }
      // mark rows
      int nMaster;
      CoinFillN(rowBlock, numberRows, -2);
      if (best2[2] < best2[0] || best2[2] < best2[1]) {
        int iRow1;
        int iRow2;
        if (best2[0] > best2[1]) {
          // Bottom rows in master
          iRow1 = row2[0] + 1;
          iRow2 = numberRows;
        } else {
          // Top rows in master
          iRow1 = 0;
          iRow2 = numberRows - row2[1];
        }
        nMaster = iRow2 - iRow1;
        CoinFillN(rowBlock + iRow1, nMaster, -1);
      } else {
        // sorted
        // Bottom rows in master (in order)
        int iRow1 = row2[2] + 1;
        nMaster = numberRows - iRow1;
        for (int i = iRow1; i < numberRows; i++)
          rowBlock[stack[i]] = -1;
      }
      if (nMaster * 2 > numberRows) {
        sprintf(generalPrint, "%d rows out of %d would be in master - no good",
          nMaster, numberRows);
        handler_->message(COIN_GENERAL_WARNING, messages_) << generalPrint << CoinMessageEol;
        delete[] rowBlock;
        delete[] columnBlock;
        delete[] whichRow;
        delete[] whichColumn;
        delete[] stack;
        CoinModel model(numberRows, numberColumns, &matrix, rowLower, rowUpper,
          columnLower, columnUpper, objective);
        model.setObjectiveOffset(objectiveOffset);
        addBlock("row_master", "column_master", model);
        return 0;
      }
    } else {
      for (iRow = 0; iRow < numberRows; iRow++)
        rowBlock[iRow] = -2;
      if (starts) {
        for (iRow = starts[0]; iRow < starts[1]; iRow++)
          rowBlock[iRow] = -1;
        for (int iBlock = 0; iBlock < maxBlocks; iBlock++) {
          for (iRow = starts[iBlock + 2]; iRow < starts[iBlock + 3]; iRow++)
            rowBlock[iRow] = iBlock;
        }
        numberBlocks = maxBlocks;
        CoinFillN(columnBlock, numberColumns, -1);
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
          CoinBigIndex kstart = columnStart[iColumn];
          CoinBigIndex kend = columnStart[iColumn] + columnLength[iColumn];
          CoinBigIndex j;
          int jBlock = -1;
          for (j = kstart; j < kend; j++) {
            int iRow = row[j];
            if (rowBlock[iRow] != -1) {
              if (jBlock == -1) {
                jBlock = rowBlock[iRow];
                columnBlock[iColumn] = jBlock;
              } else {
                assert(rowBlock[iRow] == jBlock);
              }
            }
          }
        }
        // these are master rows
      } else if (numberRows == 105127) {
        // ken-18
        for (iRow = 104976; iRow < numberRows; iRow++)
          rowBlock[iRow] = -1;
      } else if (numberRows == 2426) {
        // ken-7
        for (iRow = 2401; iRow < numberRows; iRow++)
          rowBlock[iRow] = -1;
      } else if (numberRows == 810) {
        for (iRow = 81; iRow < 84; iRow++)
          rowBlock[iRow] = -1;
      } else if (numberRows == 5418) {
        for (iRow = 564; iRow < 603; iRow++)
          rowBlock[iRow] = -1;
      } else if (numberRows == 10280) {
        // osa-60
        for (iRow = 10198; iRow < 10280; iRow++)
          rowBlock[iRow] = -1;
      } else if (numberRows == 1503) {
        // degen3
        for (iRow = 0; iRow < 561; iRow++)
          rowBlock[iRow] = -1;
      } else if (numberRows == 929) {
        // czprob
        for (iRow = 0; iRow < 39; iRow++)
          rowBlock[iRow] = -1;
      } else if (numberRows == 33874) {
        // pds-20
        for (iRow = 31427; iRow < 33874; iRow++)
          rowBlock[iRow] = -1;
      } else if (numberRows == 24902) {
        // allgrade
        int kRow = 818;
        for (iRow = 0; iRow < kRow; iRow++)
          rowBlock[iRow] = -1;
      }
    }
    if (!starts) {
      numberBlocks = 0;
      CoinFillN(columnBlock, numberColumns, -2);
    }
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      CoinBigIndex kstart = columnStart[iColumn];
      CoinBigIndex kend = columnStart[iColumn] + columnLength[iColumn];
      if (columnBlock[iColumn] == -2) {
        // column not allocated
        CoinBigIndex j;
        int nstack = 0;
        for (j = kstart; j < kend; j++) {
          int iRow = row[j];
          if (rowBlock[iRow] != -1) {
            assert(rowBlock[iRow] == -2);
            rowBlock[iRow] = numberBlocks; // mark
            stack[nstack++] = iRow;
          }
        }
        if (nstack) {
          // new block - put all connected in
          numberBlocks++;
          columnBlock[iColumn] = numberBlocks - 1;
          while (nstack) {
            int iRow = stack[--nstack];
            CoinBigIndex k;
            for (k = rowStart[iRow]; k < rowStart[iRow] + rowLength[iRow]; k++) {
              int iColumn = column[k];
              CoinBigIndex kkstart = columnStart[iColumn];
              CoinBigIndex kkend = kkstart + columnLength[iColumn];
              if (columnBlock[iColumn] == -2) {
                columnBlock[iColumn] = numberBlocks - 1; // mark
                // column not allocated
                CoinBigIndex jj;
                for (jj = kkstart; jj < kkend; jj++) {
                  int jRow = row[jj];
                  if (rowBlock[jRow] == -2) {
                    rowBlock[jRow] = numberBlocks - 1;
                    stack[nstack++] = jRow;
                  }
                }
              } else {
                assert(columnBlock[iColumn] == numberBlocks - 1);
              }
            }
          }
        } else {
          // Only in master
          columnBlock[iColumn] = -1;
        }
      }
    }
    delete[] stack;
    int numberMasterRows = 0;
    for (iRow = 0; iRow < numberRows; iRow++) {
      int iBlock = rowBlock[iRow];
      if (iBlock == -1)
        numberMasterRows++;
    }
    int numberMasterColumns = 0;
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      int iBlock = columnBlock[iColumn];
      if (iBlock == -1)
        numberMasterColumns++;
    }
    if (numberBlocks <= maxBlocks)
      sprintf(generalPrint, "%d blocks found - %d rows, %d columns in master",
        numberBlocks, numberMasterRows, numberMasterColumns);
    else
      sprintf(generalPrint, "%d blocks found (reduced to %d) - %d rows, %d columns in master",
        numberBlocks, maxBlocks, numberMasterRows, numberMasterColumns);
    handler_->message(COIN_GENERAL_INFO, messages_) << generalPrint << CoinMessageEol;
    if (numberBlocks) {
      if (numberBlocks > maxBlocks) {
        int iBlock;
        for (iRow = 0; iRow < numberRows; iRow++) {
          iBlock = rowBlock[iRow];
          if (iBlock >= 0)
            rowBlock[iRow] = iBlock % maxBlocks;
        }
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
          iBlock = columnBlock[iColumn];
          if (iBlock >= 0)
            columnBlock[iColumn] = iBlock % maxBlocks;
        }
        numberBlocks = maxBlocks;
      }
    }
    // make up problems
    // Create all sub problems
    // Space for creating
    double *obj = new double[numberColumns];
    double *columnLo = new double[numberColumns];
    double *columnUp = new double[numberColumns];
    double *rowLo = new double[numberRows];
    double *rowUp = new double[numberRows];
    // Counts
    int *rowCount = reinterpret_cast< int * >(rowLo);
    CoinZeroN(rowCount, numberBlocks);
    for (int i = 0; i < numberRows; i++) {
      int iBlock = rowBlock[i];
      if (iBlock >= 0)
        rowCount[iBlock]++;
    }
    // allocate empty rows
    for (int i = 0; i < numberRows; i++) {
      int iBlock = rowBlock[i];
      if (iBlock == -2) {
        // find block with smallest count
        int iSmall = -1;
        int smallest = numberRows;
        for (int j = 0; j < numberBlocks; j++) {
          if (rowCount[j] < smallest) {
            iSmall = j;
            smallest = rowCount[j];
          }
        }
        rowBlock[i] = iSmall;
        rowCount[iSmall]++;
      }
    }
    int *columnCount = reinterpret_cast< int * >(rowUp);
    CoinZeroN(columnCount, numberBlocks);
    for (int i = 0; i < numberColumns; i++) {
      int iBlock = columnBlock[i];
      if (iBlock >= 0)
        columnCount[iBlock]++;
    }
    int maximumSize = 0;
    for (int i = 0; i < numberBlocks; i++) {
      sprintf(generalPrint, "Block %d has %d rows and %d columns", i,
        rowCount[i], columnCount[i]);
      handler_->message(COIN_GENERAL_INFO2, messages_) << generalPrint << CoinMessageEol;
      int k = 2 * rowCount[i] + columnCount[i];
      maximumSize = std::max(maximumSize, k);
    }
    if (maximumSize * 10 > 4 * (2 * numberRows + numberColumns) && !wantDecomposition) {
      // No good
      sprintf(generalPrint, "Doesn't look good");
      handler_->message(COIN_GENERAL_WARNING, messages_) << generalPrint << CoinMessageEol;
      delete[] rowBlock;
      delete[] columnBlock;
      delete[] whichRow;
      delete[] whichColumn;
      delete[] obj;
      delete[] columnLo;
      delete[] columnUp;
      delete[] rowLo;
      delete[] rowUp;
      CoinModel model(numberRows, numberColumns, &matrix, rowLower, rowUpper,
        columnLower, columnUpper, objective);
      model.setObjectiveOffset(objectiveOffset);
      addBlock("row_master", "column_master", model);
      return 0;
    }
    // Name for master so at top
    addRowBlock(numberMasterRows, "row_master");
    // Arrays
    // get full matrix
    CoinPackedMatrix fullMatrix = matrix;
    int numberRow2, numberColumn2;
    int iBlock;
    for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
      char rowName[20];
      sprintf(rowName, "row_%d", iBlock);
      char columnName[20];
      sprintf(columnName, "column_%d", iBlock);
      numberRow2 = 0;
      numberColumn2 = 0;
      for (iRow = 0; iRow < numberRows; iRow++) {
        if (iBlock == rowBlock[iRow]) {
          rowLo[numberRow2] = rowLower[iRow];
          rowUp[numberRow2] = rowUpper[iRow];
          whichRow[numberRow2++] = iRow;
        }
      }
      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        if (iBlock == columnBlock[iColumn]) {
          obj[numberColumn2] = objective[iColumn];
          columnLo[numberColumn2] = columnLower[iColumn];
          columnUp[numberColumn2] = columnUpper[iColumn];
          whichColumn[numberColumn2++] = iColumn;
        }
      }
      // Diagonal block
      CoinPackedMatrix mat(fullMatrix,
        numberRow2, whichRow,
        numberColumn2, whichColumn);
      // make sure correct dimensions
      mat.setDimensions(numberRow2, numberColumn2);
      CoinModel *block = new CoinModel(numberRow2, numberColumn2, &mat,
        rowLo, rowUp, NULL, NULL, NULL);
      block->setOriginalIndices(whichRow, whichColumn);
      addBlock(rowName, columnName, block); // takes ownership
      // and top block
      numberRow2 = 0;
      // get top matrix
      for (iRow = 0; iRow < numberRows; iRow++) {
        int iBlock = rowBlock[iRow];
        if (iBlock == -1) {
          whichRow[numberRow2++] = iRow;
        }
      }
      CoinPackedMatrix top(fullMatrix,
        numberRow2, whichRow,
        numberColumn2, whichColumn);
      // make sure correct dimensions
      top.setDimensions(numberRow2, numberColumn2);
      block = new CoinModel(numberMasterRows, numberColumn2, &top,
        NULL, NULL, columnLo, columnUp, obj);
      block->setOriginalIndices(whichRow, whichColumn);
      addBlock("row_master", columnName, block); // takes ownership
    }
    // and master
    numberRow2 = 0;
    numberColumn2 = 0;
    for (iRow = 0; iRow < numberRows; iRow++) {
      int iBlock = rowBlock[iRow];
      if (iBlock == -1) {
        rowLo[numberRow2] = rowLower[iRow];
        rowUp[numberRow2] = rowUpper[iRow];
        whichRow[numberRow2++] = iRow;
      }
    }
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      int iBlock = columnBlock[iColumn];
      if (iBlock < 0) {
        obj[numberColumn2] = objective[iColumn];
        columnLo[numberColumn2] = columnLower[iColumn];
        columnUp[numberColumn2] = columnUpper[iColumn];
        whichColumn[numberColumn2++] = iColumn;
      }
    }
    delete[] rowBlock;
    delete[] columnBlock;
    CoinPackedMatrix top(fullMatrix,
      numberRow2, whichRow,
      numberColumn2, whichColumn);
    // make sure correct dimensions
    top.setDimensions(numberRow2, numberColumn2);
    CoinModel *block = new CoinModel(numberRow2, numberColumn2, &top,
      rowLo, rowUp,
      columnLo, columnUp, obj);
    block->setOriginalIndices(whichRow, whichColumn);
    addBlock("row_master", "column_master", block); // takes ownership
    delete[] whichRow;
    delete[] whichColumn;
    delete[] obj;
    delete[] columnLo;
    delete[] columnUp;
    delete[] rowLo;
    delete[] rowUp;
  } else if (type == 2) {
    // Try master at beginning and end
    // get row copy
    CoinPackedMatrix rowCopy = matrix;
    rowCopy.reverseOrdering();
    const int *row = matrix.getIndices();
    const int *columnLength = matrix.getVectorLengths();
    const CoinBigIndex *columnStart = matrix.getVectorStarts();
    //const double * elementByColumn = matrix.getElements();
    const int *column = rowCopy.getIndices();
    const int *rowLength = rowCopy.getVectorLengths();
    const CoinBigIndex *rowStart = rowCopy.getVectorStarts();
    //const double * elementByRow = rowCopy.getElements();
    int numberColumns = matrix.getNumCols();
    int *columnBlock = new int[numberColumns + 1];
    int iColumn;
    // Column counts (maybe look at long columns for master)
    CoinZeroN(columnBlock, numberColumns + 1);
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      int length = columnLength[iColumn];
      columnBlock[length]++;
    }
    if (!starts) {
      for (iColumn = 0; iColumn < numberColumns + 1; iColumn++) {
        if (columnBlock[iColumn]) {
          sprintf(generalPrint, "%d columns have %d elements", columnBlock[iColumn], iColumn);
          handler_->message(COIN_GENERAL_INFO2, messages_) << generalPrint << CoinMessageEol;
        }
      }
    }
    bool newWay = true;
    if (maxBlocks > 1000000) {
      newWay = false;
      maxBlocks -= 1000000;
    }
    // to say if row looked at
    int numberRows = matrix.getNumRows();
    int *rowBlock = new int[numberRows];
    int iRow;
    int *whichRow = new int[numberRows];
    int *whichColumn = new int[numberColumns];
    int *stack = new int[numberColumns];
    if (newWay) {
      //double best2[3]={COIN_DBL_MAX,COIN_DBL_MAX,COIN_DBL_MAX};
      double best2[3] = { 0.0, 0.0, 0.0 };
      int column2[3] = { -1, -1, -1 };
      // try forward and backward and sorted
      for (int iWay = 0; iWay < 3; iWay++) {
        if (iWay == 0) {
          // forwards
          for (int i = 0; i < numberColumns; i++)
            stack[i] = i;
        } else if (iWay == 1) {
          // backwards
          for (int i = 0; i < numberColumns; i++)
            stack[i] = numberColumns - 1 - i;
        } else {
          // sparsest first
          for (int i = 0; i < numberColumns; i++) {
            columnBlock[i] = columnLength[i];
            stack[i] = i;
          }
          CoinSort_2(columnBlock, columnBlock + numberColumns, stack);
        }
        CoinFillN(columnBlock, numberColumns, -1);
        columnBlock[numberColumns] = 0;
        CoinFillN(whichColumn, numberColumns, -1);
        CoinFillN(rowBlock, numberRows, -1);
        CoinFillN(whichRow, numberRows, -1);
        int numberMarkedRows = 0;
        numberBlocks = 0;
        int bestColumn = -1;
        int maximumInBlock = 0;
        int columnsDone = 0;
        int checkAfterColumns = (5 * numberColumns) / 10 + 1;
#ifdef OSL_WAY
        double best = COIN_DBL_MAX;
        int bestColumnsDone = -1;
#else
        double best = 0.0; //COIN_DBL_MAX;
#endif
        int numberGoodBlocks = 0;

        for (int kColumn = 0; kColumn < numberColumns; kColumn++) {
          iColumn = stack[kColumn];
          CoinBigIndex start = columnStart[iColumn];
          CoinBigIndex end = start + columnLength[iColumn];
          int iBlock = -1;
          for (CoinBigIndex j = start; j < end; j++) {
            int iRow = row[j];
            if (rowBlock[iRow] >= 0) {
              // already marked
              if (iBlock < 0) {
                iBlock = rowBlock[iRow];
              } else if (iBlock != rowBlock[iRow]) {
                // join two blocks
                int jBlock = rowBlock[iRow];
                numberGoodBlocks--;
                // Increase count of iBlock
                columnBlock[iBlock] += columnBlock[jBlock];
                columnBlock[jBlock] = 0;
                // First row of block jBlock
                int jRow = whichColumn[jBlock];
                while (jRow >= 0) {
                  rowBlock[jRow] = iBlock;
                  iRow = jRow;
                  jRow = whichRow[jRow];
                }
                whichRow[iRow] = whichColumn[iBlock];
                whichColumn[iBlock] = whichColumn[jBlock];
                whichColumn[jBlock] = -1;
              }
            }
          }
          int n = static_cast< int >(end - start);
          // If not in block - then start one
          if (iBlock < 0) {
            // unless null column
            if (n) {
              iBlock = numberBlocks;
              numberBlocks++;
              numberGoodBlocks++;
              int jRow = row[start];
              rowBlock[jRow] = iBlock;
              whichColumn[iBlock] = jRow;
              numberMarkedRows += n;
              columnBlock[iBlock] = n;
              for (CoinBigIndex j = start + 1; j < end; j++) {
                int iRow = row[j];
                rowBlock[iRow] = iBlock;
                whichRow[jRow] = iRow;
                jRow = iRow;
              }
            } else {
              columnBlock[numberColumns]++;
            }
          } else {
            // add all to this block if not already in
            int jRow = whichColumn[iBlock];
            for (CoinBigIndex j = start; j < end; j++) {
              int iRow = row[j];
              if (rowBlock[iRow] < 0) {
                numberMarkedRows++;
                columnBlock[iBlock]++;
                rowBlock[iRow] = iBlock;
                whichRow[iRow] = jRow;
                jRow = iRow;
              }
            }
            whichColumn[iBlock] = jRow;
          }
          columnsDone++;
          if (iBlock >= 0)
            maximumInBlock = std::max(maximumInBlock, columnBlock[iBlock]);
          if (columnsDone >= checkAfterColumns) {
            assert(numberGoodBlocks > 0);
            double averageSize = static_cast< double >(numberMarkedRows) / static_cast< double >(numberGoodBlocks);
#ifndef OSL_WAY
            double notionalBlocks = static_cast< double >(numberMarkedRows) / averageSize;
            if (maximumInBlock < 3 * averageSize && numberGoodBlocks > 2) {
              if (best * (numberColumns - columnsDone) < notionalBlocks) {
                best = notionalBlocks / static_cast< double >(numberColumns - columnsDone);
                bestColumn = kColumn;
              }
            }
#else
            if (maximumInBlock * 10 < numberRows * 11 && numberGoodBlocks > 1) {
              double test = maximumInBlock + 0.0 * averageSize;
              if (best * static_cast< double >(columnsDone) > test) {
                best = test / static_cast< double >(columnsDone);
                bestColumn = kColumn;
                bestColumnsDone = columnsDone;
              }
            }
#endif
          }
        }
#ifndef OSL_WAY
        best2[iWay] = best;
#else
        if (bestColumnsDone < numberColumns)
          best2[iWay] = -(numberColumns - bestColumnsDone);
        else
          best2[iWay] = -numberColumns;
#endif
        column2[iWay] = bestColumn;
      }
      // mark columns
      int nMaster;
      CoinFillN(columnBlock, numberColumns, -2);
      if (best2[2] < best2[0] || best2[2] < best2[1]) {
        int iColumn1;
        int iColumn2;
        if (best2[0] > best2[1]) {
          // End columns in master
          iColumn1 = column2[0] + 1;
          iColumn2 = numberColumns;
        } else {
          // Beginning columns in master
          iColumn1 = 0;
          iColumn2 = numberColumns - column2[1];
        }
        nMaster = iColumn2 - iColumn1;
        CoinFillN(columnBlock + iColumn1, nMaster, -1);
      } else {
        // sorted
        // End columns in master (in order)
        int iColumn1 = column2[2] + 1;
        nMaster = numberColumns - iColumn1;
        for (int i = iColumn1; i < numberColumns; i++)
          columnBlock[stack[i]] = -1;
      }
      if (nMaster * 2 > numberColumns) {
        sprintf(generalPrint, "%d columns out of %d would be in master - no good",
          nMaster, numberColumns);
        handler_->message(COIN_GENERAL_WARNING, messages_) << generalPrint << CoinMessageEol;
        delete[] rowBlock;
        delete[] columnBlock;
        delete[] whichRow;
        delete[] whichColumn;
        delete[] stack;
        CoinModel model(numberRows, numberColumns, &matrix, rowLower, rowUpper,
          columnLower, columnUpper, objective);
        model.setObjectiveOffset(objectiveOffset);
        addBlock("row_master", "column_master", model);
        return 0;
      }
    } else {
      for (iColumn = 0; iColumn < numberColumns; iColumn++)
        columnBlock[iColumn] = -2;
      if (starts) {
        for (iColumn = starts[0]; iColumn < starts[1]; iColumn++)
          columnBlock[iColumn] = -1;
        for (int iBlock = 0; iBlock < maxBlocks; iBlock++) {
          for (iColumn = starts[iBlock + 2]; iColumn < starts[iBlock + 3]; iColumn++)
            columnBlock[iColumn] = iBlock;
        }
        numberBlocks = maxBlocks;
        CoinFillN(rowBlock, numberRows, -1);
        for (iRow = 0; iRow < numberRows; iRow++) {
          CoinBigIndex kstart = rowStart[iRow];
          CoinBigIndex kend = rowStart[iRow] + rowLength[iRow];
          CoinBigIndex j;
          int jBlock = -1;
          for (j = kstart; j < kend; j++) {
            int iColumn = column[j];
            if (columnBlock[iColumn] != -1) {
              if (jBlock == -1) {
                jBlock = columnBlock[iColumn];
                rowBlock[iRow] = jBlock;
              } else {
                assert(columnBlock[iColumn] == jBlock);
              }
            }
          }
        }
        // these are master columns
      } else if (numberColumns == 2426) {
        // ken-7 dual
        for (iColumn = 2401; iColumn < numberColumns; iColumn++)
          columnBlock[iColumn] = -1;
      } else if (numberColumns == 10193) {
        // stormG2-8
        for (iColumn = 0; iColumn < 121; iColumn++)
          columnBlock[iColumn] = -1;
      } else if (numberColumns == 157496) {
        // stormG2-125
        for (iColumn = 0; iColumn < 121; iColumn++)
          columnBlock[iColumn] = -1;
      } else if (numberColumns == 1259121) {
        // stormG2_1000
        for (iColumn = 0; iColumn < 121; iColumn++)
          columnBlock[iColumn] = -1;
      }
    }
    if (!starts) {
      numberBlocks = 0;
      CoinFillN(rowBlock, numberRows, -2);
    }
    for (iRow = 0; iRow < numberRows; iRow++) {
      CoinBigIndex kstart = rowStart[iRow];
      CoinBigIndex kend = rowStart[iRow] + rowLength[iRow];
      if (rowBlock[iRow] == -2) {
        // row not allocated
        CoinBigIndex j;
        int nstack = 0;
        for (j = kstart; j < kend; j++) {
          int iColumn = column[j];
          if (columnBlock[iColumn] != -1) {
            assert(columnBlock[iColumn] == -2);
            columnBlock[iColumn] = numberBlocks; // mark
            stack[nstack++] = iColumn;
          }
        }
        if (nstack) {
          // new block - put all connected in
          numberBlocks++;
          rowBlock[iRow] = numberBlocks - 1;
          while (nstack) {
            int iColumn = stack[--nstack];
            CoinBigIndex k;
            for (k = columnStart[iColumn]; k < columnStart[iColumn] + columnLength[iColumn]; k++) {
              int iRow = row[k];
              CoinBigIndex kkstart = rowStart[iRow];
              CoinBigIndex kkend = kkstart + rowLength[iRow];
              if (rowBlock[iRow] == -2) {
                rowBlock[iRow] = numberBlocks - 1; // mark
                // row not allocated
                CoinBigIndex jj;
                for (jj = kkstart; jj < kkend; jj++) {
                  int jColumn = column[jj];
                  if (columnBlock[jColumn] == -2) {
                    columnBlock[jColumn] = numberBlocks - 1;
                    stack[nstack++] = jColumn;
                  }
                }
              } else {
                assert(rowBlock[iRow] == numberBlocks - 1);
              }
            }
          }
        } else {
          // Only in master
          rowBlock[iRow] = -1;
        }
      }
    }
    delete[] stack;
    int numberMasterColumns = 0;
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      int iBlock = columnBlock[iColumn];
      if (iBlock == -1)
        numberMasterColumns++;
    }
    int numberMasterRows = 0;
    for (iRow = 0; iRow < numberRows; iRow++) {
      int iBlock = rowBlock[iRow];
      if (iBlock == -1)
        numberMasterRows++;
    }
    if (numberBlocks <= maxBlocks)
      sprintf(generalPrint, "%d blocks found - %d columns, %d rows in master",
        numberBlocks, numberMasterColumns, numberMasterRows);
    else
      sprintf(generalPrint, "%d blocks found (reduced to %d) - %d columns, %d rows in master",
        numberBlocks, maxBlocks, numberMasterColumns, numberMasterRows);
    handler_->message(COIN_GENERAL_INFO, messages_) << generalPrint << CoinMessageEol;
    if (numberBlocks) {
      if (numberBlocks > maxBlocks) {
        int iBlock;
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
          iBlock = columnBlock[iColumn];
          if (iBlock >= 0)
            columnBlock[iColumn] = iBlock % maxBlocks;
        }
        for (iRow = 0; iRow < numberRows; iRow++) {
          iBlock = rowBlock[iRow];
          if (iBlock >= 0)
            rowBlock[iRow] = iBlock % maxBlocks;
        }
        numberBlocks = maxBlocks;
      }
    }
    // make up problems
    // Create all sub problems
    // Space for creating
    double *obj = new double[numberColumns];
    double *rowLo = new double[numberRows];
    double *rowUp = new double[numberRows];
    double *columnLo = new double[numberColumns];
    double *columnUp = new double[numberColumns];
    // Counts
    int *columnCount = reinterpret_cast< int * >(columnLo);
    CoinZeroN(columnCount, numberBlocks);
    for (int i = 0; i < numberColumns; i++) {
      int iBlock = columnBlock[i];
      if (iBlock >= 0)
        columnCount[iBlock]++;
    }
    // allocate empty columns
    for (int i = 0; i < numberColumns; i++) {
      int iBlock = columnBlock[i];
      if (iBlock == -2) {
        // find block with smallest count
        int iSmall = -1;
        int smallest = numberColumns;
        for (int j = 0; j < numberBlocks; j++) {
          if (columnCount[j] < smallest) {
            iSmall = j;
            smallest = columnCount[j];
          }
        }
        columnBlock[i] = iSmall;
        columnCount[iSmall]++;
      }
    }
    int *rowCount = reinterpret_cast< int * >(columnUp);
    CoinZeroN(rowCount, numberBlocks);
    for (int i = 0; i < numberRows; i++) {
      int iBlock = rowBlock[i];
      if (iBlock >= 0)
        rowCount[iBlock]++;
    }
    int maximumSize = 0;
    for (int i = 0; i < numberBlocks; i++) {
      sprintf(generalPrint, "Block %d has %d columns and %d rows", i,
        columnCount[i], rowCount[i]);
      handler_->message(COIN_GENERAL_INFO2, messages_) << generalPrint << CoinMessageEol;
      int k = 2 * columnCount[i] + rowCount[i];
      maximumSize = std::max(maximumSize, k);
    }
    if ((maximumSize * 10 > 4 * (2 * numberColumns + numberRows) || numberMasterRows * 10 > numberRows) && !wantDecomposition) {
      // No good
      sprintf(generalPrint, "Doesn't look good");
      handler_->message(COIN_GENERAL_WARNING, messages_) << generalPrint << CoinMessageEol;
      delete[] rowBlock;
      delete[] columnBlock;
      delete[] whichRow;
      delete[] whichColumn;
      delete[] obj;
      delete[] columnLo;
      delete[] columnUp;
      delete[] rowLo;
      delete[] rowUp;
#if 0
      CoinModel model(numberRows,numberColumns,&matrix, rowLower, rowUpper,
		      columnLower,columnUpper,objective);
      model.setObjectiveOffset(objectiveOffset);
      addBlock("row_master","column_master",model);
#endif
      return 0;
    }
    // Name for master so at beginning
    addColumnBlock(numberMasterColumns, "column_master");
    // Arrays
    // get full matrix
    CoinPackedMatrix fullMatrix = matrix;
    int numberRow2, numberColumn2;
    int iBlock;
    for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
      char rowName[20];
      sprintf(rowName, "row_%d", iBlock);
      char columnName[20];
      sprintf(columnName, "column_%d", iBlock);
      numberRow2 = 0;
      numberColumn2 = 0;
      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        if (iBlock == columnBlock[iColumn]) {
          obj[numberColumn2] = objective[iColumn];
          columnLo[numberColumn2] = columnLower[iColumn];
          columnUp[numberColumn2] = columnUpper[iColumn];
          whichColumn[numberColumn2++] = iColumn;
        }
      }
      for (iRow = 0; iRow < numberRows; iRow++) {
        if (iBlock == rowBlock[iRow]) {
          rowLo[numberRow2] = rowLower[iRow];
          rowUp[numberRow2] = rowUpper[iRow];
          whichRow[numberRow2++] = iRow;
        }
      }
      // Diagonal block
      CoinPackedMatrix mat(fullMatrix,
        numberRow2, whichRow,
        numberColumn2, whichColumn);
      // make sure correct dimensions
      mat.setDimensions(numberRow2, numberColumn2);
      CoinModel *block = new CoinModel(numberRow2, numberColumn2, &mat,
        rowLo, rowUp, columnLo, columnUp, obj);
      block->setOriginalIndices(whichRow, whichColumn);
      addBlock(rowName, columnName, block); // takes ownership
      // and beginning block
      numberColumn2 = 0;
      // get beginning matrix
      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        int iBlock = columnBlock[iColumn];
        if (iBlock == -1) {
          whichColumn[numberColumn2++] = iColumn;
        }
      }
      CoinPackedMatrix beginning(fullMatrix,
        numberRow2, whichRow,
        numberColumn2, whichColumn);
      // make sure correct dimensions *********
      beginning.setDimensions(numberRow2, numberColumn2);
      block = new CoinModel(numberRow2, numberMasterColumns, &beginning,
        NULL, NULL, NULL, NULL, NULL);
      block->setOriginalIndices(whichRow, whichColumn);
      addBlock(rowName, "column_master", block); // takes ownership
    }
    // and master
    numberRow2 = 0;
    numberColumn2 = 0;
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
      int iBlock = columnBlock[iColumn];
      if (iBlock == -1) {
        obj[numberColumn2] = objective[iColumn];
        columnLo[numberColumn2] = columnLower[iColumn];
        columnUp[numberColumn2] = columnUpper[iColumn];
        whichColumn[numberColumn2++] = iColumn;
      }
    }
    for (iRow = 0; iRow < numberRows; iRow++) {
      int iBlock = rowBlock[iRow];
      if (iBlock < 0) {
        rowLo[numberRow2] = rowLower[iRow];
        rowUp[numberRow2] = rowUpper[iRow];
        whichRow[numberRow2++] = iRow;
      }
    }
    delete[] rowBlock;
    delete[] columnBlock;
    CoinPackedMatrix beginning(fullMatrix,
      numberRow2, whichRow,
      numberColumn2, whichColumn);
    // make sure correct dimensions
    beginning.setDimensions(numberRow2, numberColumn2);
    CoinModel *block = new CoinModel(numberRow2, numberColumn2, &beginning,
      rowLo, rowUp,
      columnLo, columnUp, obj);
    block->setOriginalIndices(whichRow, whichColumn);
    addBlock("row_master", "column_master", block); // takes ownership
    delete[] whichRow;
    delete[] whichColumn;
    delete[] obj;
    delete[] columnLo;
    delete[] columnUp;
    delete[] rowLo;
    delete[] rowUp;
  } else {
    abort();
  }
  return numberBlocks;
}
/* Decompose a CoinModel
   1 - try D-W
   2 - try Benders
   3 - try Staircase
   Returns number of blocks or zero if no structure
*/
int CoinStructuredModel::decompose(const CoinModel &coinModel, int type,
  int maxBlocks, const char **starts)
{
  const CoinPackedMatrix *matrix = coinModel.packedMatrix();
  assert(matrix != NULL);
  // Arrays
  const double *objective = coinModel.objectiveArray();
  const double *columnLower = coinModel.columnLowerArray();
  const double *columnUpper = coinModel.columnUpperArray();
  const double *rowLower = coinModel.rowLowerArray();
  const double *rowUpper = coinModel.rowUpperArray();
  int *blockStarts = NULL;
  char generalPrint[200];
  bool wantDecomposition = maxBlocks > 1; // flag to say we really want it decomposed
  int numberTotal = coinModel.numberColumns() + coinModel.numberRows();
  if (maxBlocks < 2 || 2 * maxBlocks > numberTotal) {
    // allow at least 400 per problem
    maxBlocks = (numberTotal + 399) / 400;
    // but gate at 8 and 1000
    maxBlocks = std::max(8, std::min(maxBlocks, 1000));
    sprintf(generalPrint, "Trying for %d blocks", maxBlocks);
    handler_->message(COIN_GENERAL_INFO, messages_) << generalPrint << CoinMessageEol;
  }
  if (starts) {
    assert(type < 3);
    // maxBlocks is number of blocks
    // first two are master - then starts
    // master at beginning or end
    blockStarts = new int[maxBlocks + 3];
    if (type == 2) {
      // column names
      int numberColumns = coinModel.numberColumns();
      // Do master block
      int iColumn;
      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        if (!strcmp(starts[0], coinModel.getColumnName(iColumn)))
          break;
      }
      if (iColumn == numberColumns) {
        sprintf(generalPrint, "Unable to find start of master block %s", starts[0]);
        handler_->message(COIN_GENERAL_WARNING, messages_) << generalPrint << CoinMessageEol;
        delete[] blockStarts;
        return 0;
      }
      if (iColumn == 0) {
        blockStarts[0] = 0;
        blockStarts[1] = -1;
      } else {
        blockStarts[0] = iColumn;
        blockStarts[1] = numberColumns;
      }
      int nBlocks = 2;
      for (int iBlock = 1; iBlock < maxBlocks + 1; iBlock++) {
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
          if (!strcmp(starts[iBlock], coinModel.getColumnName(iColumn)))
            break;
        }
        if (iColumn == numberColumns) {
          sprintf(generalPrint, "Unable to find start of block %d %s", iBlock, starts[iBlock]);
          handler_->message(COIN_GENERAL_WARNING, messages_) << generalPrint << CoinMessageEol;
          delete[] blockStarts;
          return 0;
        }
        if (nBlocks == 2 && blockStarts[1] == -1)
          blockStarts[1] = iColumn - 1;
        blockStarts[nBlocks++] = iColumn;
      }
      if (blockStarts[1] == numberColumns)
        blockStarts[nBlocks++] = blockStarts[0];
      else
        blockStarts[nBlocks++] = numberColumns;
    } else {
      // row names
      int numberRows = coinModel.numberRows();
      // Do master block
      int iRow;
      for (iRow = 0; iRow < numberRows; iRow++) {
        if (!strcmp(starts[0], coinModel.getRowName(iRow)))
          break;
      }
      if (iRow == numberRows) {
        sprintf(generalPrint, "Unable to find start of master block %s", starts[0]);
        handler_->message(COIN_GENERAL_WARNING, messages_) << generalPrint << CoinMessageEol;
        delete[] blockStarts;
        return 0;
      }
      if (iRow == 0) {
        blockStarts[0] = 0;
        blockStarts[1] = -1;
      } else {
        blockStarts[0] = iRow;
        blockStarts[1] = numberRows;
      }
      int nBlocks = 2;
      for (int iBlock = 1; iBlock < maxBlocks + 1; iBlock++) {
        for (iRow = 0; iRow < numberRows; iRow++) {
          if (!strcmp(starts[iBlock], coinModel.getRowName(iRow)))
            break;
        }
        if (iRow == numberRows) {
          sprintf(generalPrint, "Unable to find start of block %d %s", iBlock, starts[iBlock]);
          handler_->message(COIN_GENERAL_WARNING, messages_) << generalPrint << CoinMessageEol;
          delete[] blockStarts;
          return 0;
        }
        if (nBlocks == 2 && blockStarts[1] == -1)
          blockStarts[1] = iRow - 1;
        blockStarts[nBlocks++] = iRow;
      }
      if (blockStarts[1] == numberRows)
        blockStarts[nBlocks++] = blockStarts[0];
      else
        blockStarts[nBlocks++] = numberRows;
    }
  }
  if (wantDecomposition)
    type += 10;
  int nBlocks = decompose(*matrix,
    rowLower, rowUpper,
    columnLower, columnUpper,
    objective, type, maxBlocks, blockStarts,
    coinModel.objectiveOffset());
  delete[] blockStarts;
  return nBlocks;
}
// Read SMPS model
int CoinStructuredModel::readSmps(const char *filename,
  bool keepNames, bool ignoreErrors)
{
  abort();
  return 0;
}
// Return block corresponding to row and column
const CoinBaseModel *
CoinStructuredModel::block(int row, int column) const
{
  const CoinBaseModel *block = NULL;
  if (blockType_) {
    for (int iBlock = 0; iBlock < numberElementBlocks_; iBlock++) {
      if (blockType_[iBlock].rowBlock == row && blockType_[iBlock].columnBlock == column) {
        block = blocks_[iBlock];
        break;
      }
    }
  }
  return block;
}
// Return block corresponding to row and column as CoinModel
const CoinBaseModel *
CoinStructuredModel::coinBlock(int row, int column) const
{
  const CoinModel *block = NULL;
  if (blockType_) {
    for (int iBlock = 0; iBlock < numberElementBlocks_; iBlock++) {
      if (blockType_[iBlock].rowBlock == row && blockType_[iBlock].columnBlock == column) {
        block = dynamic_cast< CoinModel * >(blocks_[iBlock]);
        assert(block);
        break;
      }
    }
  }
  return block;
}
/* Fill pointers corresponding to row and column.
   False if any missing */
CoinModelBlockInfo
CoinStructuredModel::block(int row, int column,
  const double *&rowLower, const double *&rowUpper,
  const double *&columnLower, const double *&columnUpper,
  const double *&objective) const
{
  CoinModelBlockInfo info;
  //memset(&info,0,sizeof(info));
  rowLower = NULL;
  rowUpper = NULL;
  columnLower = NULL;
  columnUpper = NULL;
  objective = NULL;
  if (blockType_) {
    for (int iBlock = 0; iBlock < numberElementBlocks_; iBlock++) {
      CoinModel *thisBlock = coinBlock(iBlock);
      if (blockType_[iBlock].rowBlock == row) {
        if (blockType_[iBlock].rhs) {
          info.rhs = 1;
          rowLower = thisBlock->rowLowerArray();
          rowUpper = thisBlock->rowUpperArray();
        }
      }
      if (blockType_[iBlock].columnBlock == column) {
        if (blockType_[iBlock].bounds) {
          info.bounds = 1;
          columnLower = thisBlock->columnLowerArray();
          columnUpper = thisBlock->columnUpperArray();
          objective = thisBlock->objectiveArray();
        }
      }
    }
  }
  return info;
}
// Return block number corresponding to row and column
int CoinStructuredModel::blockIndex(int row, int column) const
{
  int block = -1;
  if (blockType_) {
    for (int iBlock = 0; iBlock < numberElementBlocks_; iBlock++) {
      if (blockType_[iBlock].rowBlock == row && blockType_[iBlock].columnBlock == column) {
        block = iBlock;
        break;
      }
    }
  }
  return block;
}

/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/
